#ifndef CONVERTER_COLOR_DIFF_CALC_H
#define CONVERTER_COLOR_DIFF_CALC_H

#include <common/color.h>
#include <cmath>

class IColorDiffCalc
{
public:
    inline double GetColorDiff(Color c1, Color c2)
    {
        auto cd = GetRGBDiff(c1.r, c1.g, c1.b, c2.r, c2.g, c2.b);

        // 此处当a较小时弱化rgb距离的结果, 当a较大时强化rgb距离结果, 同时强调a的差距
        // TODO: 效果等待验证
        // 失败, 如何计算alpha
//        return cd * (c1.a + c2.a) / 128 + 3 * (c2.a - c1.a);
        return cd;
    }

private:
    virtual double GetRGBDiff(uint8_t r1, uint8_t g1, uint8_t b1, uint8_t r2, uint8_t g2, uint8_t b2) = 0;
};

class HSVColorDiffCalc : public IColorDiffCalc
{
private:
    /**
     * @brief https://blog.csdn.net/hai29785/article/details/55548798
     */
    void RGB2HSV(uint8_t r, uint8_t g, uint8_t b, double& h, double& s, double& v)
    {
        double rD = double(r) / 255.0;
        double gD = double(g) / 255.0;
        double bD = double(b) / 255.0;
        double cMax = std::max(std::max(rD, gD), bD);
        double cMin = std::min(std::min(rD, gD), bD);
        double delta = cMax - cMin;

        if (cMax == 0)
            h = 0;
        else if (cMax == rD)
            h = 60 * (((gD - bD) / delta) + 0);
        else if (cMax == gD)
            h = 60 * (((bD - rD) / delta) + 2);
        else if (cMax == bD)
            h = 60 * (((rD - gD) / delta) + 4);

        if (cMax == 0)
            s = 0;
        else
            s = delta / cMax;

        v = cMax;
    }

    double GetRGBDiff(uint8_t r1, uint8_t g1, uint8_t b1, uint8_t r2, uint8_t g2, uint8_t b2) override
    {
        double h1, s1, v1, h2, s2, v2;
        RGB2HSV(r1, g1, b1, h1, s1, v1);
        RGB2HSV(r2, g2, b2, h2, s2, v2);
        double x1 = r1 * v1 * s1 * cos(h1 / 180 * 3.14159265358979323846264338327950288);
        double y1 = r1 * v1 * s1 * sin(h1 / 180 * 3.14159265358979323846264338327950288);
        double z1 = h1 * (1 - v1);
        double x2 = r2 * v2 * s2 * cos(h2 / 180 * 3.14159265358979323846264338327950288);
        double y2 = r2 * v2 * s2 * sin(h2 / 180 * 3.14159265358979323846264338327950288);
        double z2 = h2 * (1 - v2);
        double dx = x1 - x2;
        double dy = y1 - y2;
        double dz = z1 - z2;
        return dx * dx + dy * dy + dz * dz;
    }
};

class EulerColorDiffCalc : public IColorDiffCalc
{
private:
    double GetRGBDiff(uint8_t r1, uint8_t g1, uint8_t b1, uint8_t r2, uint8_t g2, uint8_t b2) override
    {
        auto r_ = (double)(r1 + r2) / 2.0 ;
        auto rd = (double)(r1 - r2);
        auto gd = (double)(g1 - g2);
        auto bd = (double)(b1 - b2);
        auto d = (2.0 + r_ / 256.0) * rd * rd + 4.0 * gd * gd + (2.0 + (255.0 - r_) / 256.0) * bd * bd;
        return d;
    }
};

class CIEDE2000ColorDiffCalc : public IColorDiffCalc
{
private:
    /**
     * @brief https://blog.csdn.net/lz0499/article/details/77345166
     */
    const inline static double param_16116 = 16.0f / 116.0f;
    const inline static double xn = 0.950456f;
    const inline static double yn = 1.0f;
    const inline static double zn = 1.088754f;
    static void RGB2CIELab(uint8_t red, uint8_t green, uint8_t blue, double& l, double& a, double& b)
    {
        auto gamma = [](double x)
        {
            return x > 0.04045 ? std::pow((x + 0.055f) / 1.055f,2.4f) : (x / 12.92);
        };

        double rg = gamma(red / 255.0);
        double gg = gamma(green / 255.0);
        double bg = gamma(blue / 255.0);

        auto f = [](double x)
        {
            return x > 0.008856f ? std::pow(x, 1.0f / 3.0f) : 7.787f * x + (16.0f / 116.0f);
        };

        double x = f((0.4124564 * rg + 0.3575761 * gg + 0.1804375 * bg) / xn);
        double y = f((0.2126729 * rg + 0.7151522 * gg + 0.0721750 * bg) / yn);
        double z = f((0.0193339 * rg + 0.1191920 * gg + 0.9503041 * bg) / zn);

        l = 116.0f * y - 16.0f;
        l = l > 0.0f ? l : 0.0f;
        a = 500.0f * (x - y);
        b = 200.0f * (y - z);
    }
    double GetRGBDiff(uint8_t red1, uint8_t green1, uint8_t blue1, uint8_t red2, uint8_t green2, uint8_t blue2) override
    {
        double l1, a1, b1, l2, a2, b2;
        RGB2CIELab(red1, green1, blue1, l1, a1, b1);
        RGB2CIELab(red2, green2, blue2, l2, a2, b2);

        auto deg2Rad = [](double x){ return (x * (3.14159265358979323846264338327950288 / 180.0)); };
        auto rad2Deg = [](double x){ return ((180.0 / 3.14159265358979323846264338327950288) * x); };

        const double k_L = 1.0, k_C = 1.0, k_H = 1.0;
        const double deg360InRad = deg2Rad(360.0);
        const double deg180InRad = deg2Rad(180.0);
        const double pow25To7 = 6103515625.0; /* pow(25, 7) */

        double C1 = sqrt((a1 * a1) + (b1 * b1));
        double C2 = sqrt((a2 * a2) + (b2 * b2));
        double barC = (C1 + C2) / 2.0;
        double G = 0.5 * (1 - sqrt(pow(barC, 7) / (pow(barC, 7) + pow25To7)));
        double a1Prime = (1.0 + G) * a1;
        double a2Prime = (1.0 + G) * a2;
        double CPrime1 = sqrt((a1Prime * a1Prime) + (b1 * b1));
        double CPrime2 = sqrt((a2Prime * a2Prime) + (b2 * b2));
        double hPrime1;
        if (b1 == 0 && a1Prime == 0)
        {
            hPrime1 = 0.0;
        }
        else
        {
            hPrime1 = atan2(b1, a1Prime);
            if (hPrime1 < 0)
                hPrime1 += deg360InRad;
        }
        double hPrime2;
        if (b2 == 0 && a2Prime == 0)
        {
            hPrime2 = 0.0;
        }
        else
        {
            hPrime2 = atan2(b2, a2Prime);
            if (hPrime2 < 0)
                hPrime2 += deg360InRad;
        }

        double deltaLPrime = l2 - l1;
        double deltaCPrime = CPrime2 - CPrime1;
        double deltahPrime;
        double CPrimeProduct = CPrime1 * CPrime2;
        if (CPrimeProduct == 0)
        {
            deltahPrime = 0;
        }
        else
        {
            deltahPrime = hPrime2 - hPrime1;
            if (deltahPrime < -deg180InRad)
                deltahPrime += deg360InRad;
            else if (deltahPrime > deg180InRad)
                deltahPrime -= deg360InRad;
        }
        double deltaHPrime = 2.0 * sqrt(CPrimeProduct) * sin(deltahPrime / 2.0);
        double barLPrime = (l1 + l2) / 2.0;
        double barCPrime = (CPrime1 + CPrime2) / 2.0;
        double barhPrime, hPrimeSum = hPrime1 + hPrime2;
        if (CPrime1 * CPrime2 == 0)
        {
            barhPrime = hPrimeSum;
        }
        else
        {
            if (fabs(hPrime1 - hPrime2) <= deg180InRad)
            {
                barhPrime = hPrimeSum / 2.0;
            }
            else
            {
                if (hPrimeSum < deg360InRad)
                    barhPrime = (hPrimeSum + deg360InRad) / 2.0;
                else
                    barhPrime = (hPrimeSum - deg360InRad) / 2.0;
            }
        }
        double T = 1.0 - (0.17 * cos(barhPrime - deg2Rad(30.0))) +
                   (0.24 * cos(2.0 * barhPrime)) +
                   (0.32 * cos((3.0 * barhPrime) + deg2Rad(6.0))) -
                   (0.20 * cos((4.0 * barhPrime) - deg2Rad(63.0)));
        double deltaTheta = deg2Rad(30.0) *
                            exp(-pow((barhPrime - deg2Rad(275.0)) / deg2Rad(25.0), 2.0));
        double R_C = 2.0 * sqrt(pow(barCPrime, 7.0) /
                                (pow(barCPrime, 7.0) + pow25To7));
        double S_L = 1 + ((0.015 * pow(barLPrime - 50.0, 2.0)) /
                          sqrt(20 + pow(barLPrime - 50.0, 2.0)));
        double S_C = 1 + (0.045 * barCPrime);
        double S_H = 1 + (0.015 * barCPrime * T);
        double R_T = (-sin(2.0 * deltaTheta)) * R_C;
        double deltaE =
                pow(deltaLPrime / (k_L * S_L), 2.0) +
                pow(deltaCPrime / (k_C * S_C), 2.0) +
                pow(deltaHPrime / (k_H * S_H), 2.0) +
                (R_T * (deltaCPrime / (k_C * S_C)) * (deltaHPrime / (k_H * S_H)));
        return (deltaE);
    }
};

#endif //CONVERTER_COLOR_DIFF_CALC_H
